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ABSTRACT 

Ho 1 s Algorithm is reviewed and demonstrated with ana- 
lytic examples. A digital computer program is developed to 
implement the algorithm for single— input, single-output 
systems and used to identify linear continuous and station- 
ary systems which are driven with a unit step as the test 
input. Discrete realization of the continuous systems is 
obtained using the measured output-samples, to a step in- 
put, directly in the algorithm. 
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I. 



REVIEW CF Ho 1 s ALGORITHM 



A. INTRODUCTION 

Design of a control system which will perform accept- 
ably with respect to some criterion and satisfy the pos- 
sible constraints over its operating range depends on 
knowledge of the process dynamics. The dynamics might be 
linear, nonlinear, time variant, stationary or might be a 
function of environment. In any case the dynamics of the 
process must be formulated by a set of differential equa- 
tions (if it is a continuously operating system) or by a 
set of difference equations (if it is a discrete time sys- 
tem) . An effective design and/or analysis can proceed 
after this formulation. 

These equations might be completely known or might be 
written down from the parameter values supplied by the 
manufacturer using the laws of physics. Sometimes only 
partial information or no information is supplied about 
the process; yet somehow the dynamics must be formulated. 
This leads to the solution of an identification problem. 

B. BACKGROUND 

A number of algorithms can be found Cl, 2,3,5]'*' for 
the solution of the identification problem. These range 



end. 



1 

Numbers in brackets 



indicate the references at the 
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from a partial identification, like the damping factor of a 
dominant complex pair of poles, to a complete identifica- 
tion of the systems dynamics using frequency and time do- 
main techniques. 

The identification problem can be defined in general 
as: to find an internal description of the system (a set 

of differential or difference equations) from the given 
external description (input-output relation) . 

Input-output relations may be defined in the time do- 
main (impulse response for continuous systems and pulse 
response for discrete systems) or in the frequency domain 
(transfer function). In the identification problem this 
relation is given in general as experimental data. 

Most of the algorithms are applicable to linear, sta- 
tionary systems only; and only this class of systems will 
be considered in this study. 

The following definitions and notations will be used. 

A continuous, linear, stationary system will be re- 
presented by 

x = Ax + Bu ~ 

•-'w /V / 

( 1 ) 

y = Cx 
where, 

A is a nxn system matrix 

B is a nxm distribution matrix 



2 

x denotes 



d 

dt 



x 
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c 


is a pxn 


observation matrix 


X 


= x ( t) is 


a 


nxl 


state vector 


u 

/%/ 


= u(t) is 


a 


mxl 


input vector 


y 


= y ( t) is 


a 


pxl 


output vector 



and a discrete, linear, stationary system will be repre- 
sented by 

x(k+l) = A x(k) + B u(k) 

U U (2 

y(k) = C D x(k) 



with dimensions same as (1) 

A continuous system is said to be controllable [ 2 ] if 
the matrix 



f B ! AB ! a 2 b ! 

I— • I I 



! A n 1 b] 



(3) 



has rank n; and is said to be observable if the matrix 



fc T ! a t c t • <A T ) 2 c t ! .... i (a t ) n-1 c T l 

““ I 1 .1 -J 



(4) 



has rank n. Same definitions hold for a discrete system by 
replacing A, B, C with A^, B^, respectively. 

A system may always be partitioned into four possible 
subsystems [2,3] as shown in Fig. 1. 

Part I: Controllable and observable 

Part II: Uncontrollable but observable 

Part III: Controllable but unobservable 

Part IV: Uncontrollable and unobservable. 



3 

Vertical dotted lines indicate partitioning. 

4 T 

C indicates transpose of C. 
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u 





Fig. 1 

As is seen from Fig. 1 only Part I gives the input- 
output relation. Although Part II seems to be contribut- 
ing to the output, if no energy is stored at this part 
(zero initial condition) any contribution will be the re- 
sult of a noise input only. 

If the system is of the order n and the subsystems 



are of order n^, n^, n m' n iv res P ect ively then 



n n j + n ZI + n m + n jv 



(5) 



Identification (alternatively called realization) of 
the system in Fig. 1 will give a system of order n . 

The above statements are given as the following 
theorem in [2]. "Knowledge of the impulse response identi- 
fies the completely controllable and completely observable 
part, and this part alone, of the dynamical system which 
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generated it. This part is itself a dynamical system and 
has the smallest dimension among all realizations. More- 
over, this part is identified by its impulse response 
uniquely up to algebraic equivalence". 

Therefore complete identification of a dynamical sys- 
tem necessitates complete controllability and observability 
of the system. Part I must constitute all the system and 
parts II, III and IV must be missing. 

In this study only completely controllable and observ- 
able systems will be considered. 



C. Ho 1 s ALGORITHM [l,4] 

First the algorithm will be given for a discrete sys- 
tem as represented in (2) . 

Let h(k)^j be the present response at output i to a 
unit pulse at input j applied k time units ago. Let H^ be 
the matrix 



= 



hfk) 1;L 


h(k) 12 


h(k) 21 


h (k) 22 



h (k) 



im 



h (k) 



Pi 



h (k) 



pm 



(6) 
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Let HO^ be the block matrix 



H0 1 = 



H 1 H 2 



H 2 H 3 



H, 



H, 



.H 



21-1 



(7) 



Let be the rank of HO^. Build the block matrices 



HO, 



1 = 1,2,3, ... and each time find the rank p., until 



p^ equal p^ + ^, let this value of p^ be n and SL be q; then 
find the elementary transformation matrices P and Q such 
that 



P(HO g )Q = 



1 0 . 



0 1 



n 1 ' s 



0 



,0 



( 8 ) 
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Save P and Q matrices. Let HOT be the block matrix 

q 



HOT = 

q 



H, 



H, 



H 



H- 



H, 



q+1 • 



H 



q+i 



H 



2q 



( 9 ) 



Then a state variable description of the system is 
found by letting 

A_ = the nxn northwest corner of P(HOT )Q 

D q 

B q = the nxm northwest corner of P (HO^) (10) 

C D = the pxn northwest corner of (HO^)Q 
The following results can be stated 

1. 

realization. 

2. By choosing P and Q suitably all possible realiza- 
tions may be obtained. 

3. Any two minimal realizations are isomorphic, that 
is there exists a nonsingular matrix W such that 



p = p ,. = n is the minimum dimension for the 

q q+i 



A_ = W A_ W 
D 2 D 1 



-1 



% - w % 



( 11 ) 



c = c w 

D 2 D 1 



-1 



11 



(12) 



and W is given explicitly as 

« - < V 2 V 2 T > _1 V 2 V 1 T 



where V r r = 1,2 is given by 

v r ■ [cl \ »; c T : ) 2 c I 



t . q-i t 



:< a d > 

, r 



T I 
D J 



(13) 



r _ r ! 

For continuous systems the algorithm follows the same 
pattern except in (6) each h(k)^j must be replaced by the 
impulse response and its successive derivatives evaluated 
at t = 0, that is, if g(t)^j is defined as the impulse 
response between input j and output i then 



h (1) . . = g(0) . . 
ij 3 ' ' .1 j 

h(2).. = g(0).. (14) 

h (3) . . = g(0) . . 
l j 3 l j 



or equally g ( t) ^ ^ must be expanded in a power series for 



all t as 



g ( t) i j = I a k t k V(k-l) i (15) 

k=l 

and each h(k) . . must be replaced by a, . These necessitate 
1J K 

that 9(t)^j must be a real analytic function. 

If the input-output relation is given in the trans- 
form domain (H(z)^ for discrete systems and H(s)^j for 
continuous systems) it must be representable as a power 
series in the form of 



12 



(s for cont. systems) 



(16) 



H 



(z) ii * 1 



\ z 



-k 



k=l 



and each h(k) i j in (6) must be replaced by b^, after that 
the algorithm is developed in similar fashion. 

In any case H^, , k = 1,2, ... constitutes a sequence of 

constant matrices. In a broad sense the identification 
problem may be restated as Cl]: Given a sequence of pxm 

constant matrices, H^, k = 1,2, ... find the triple 
Ap, Bjy Cp (or A, B, C) of constant matrices such that 

= C D A D k_1 B D , k = 1,2, ... (17) 

to show the mechanics of the algorithm the following ex- 
amples are presented. 

1. Example 1 Single -Input ^ Single-Output 

Discrete System (m = 1, p = 1) 



Le t the 



k = 1,2, 



be 



H 1 = 1 



H 4 = 6 



H 2 = 


0 


H 3 = 


-2 


H 5 = 


-14 


H 6 = 


30 



Start building HO^, H = 1,2, ... 

h °i - C h i ] - [ 1 j ’ Pi = 1 





" H ! 


H 2 




1 


0 


ho 2 = 






= 








H 2 


H 3_ 




0 


-2 



; P 2 2 
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HO, 



10-2 
0-2 6 
-2 6 -14 



P 3 = 2 



therefore 

n = 2 and q = 2 



HOT, 



H 2 


H 3 




0 


-2 


_ H 3 


H 4_ 




-2 


6 



A set of P and Q matrices which satisfy (8) is 



0 



/ Qi 



p 1 (hot 2 )q 1 



p 1 (ho 2 ) = 



0 



-2 -3 



0 



0 -2 



= K 






B 



D n 



(ho 2 )q 1 



0 



0 



=> C D, - [ ] 
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therefore the state equations are 



0 



x(k+l) = 



xOO + 



u(k) 



-2 -3 



0 



yOO 



- [ 1 o] S (k) 



this realization must satisfy (17) 



»k - ■ * - 1 - 2 - 



k— 1 

A can be calculated by the Cayley-Hamil ton technique 
D 1 



A* 1 = a i + a A 

Dj o 1 Dj 



ot 



ol 






-t »] 





- 2tt l 


ot 

o 


-3tt 1 


ot 

o 


“l 




1 


! 

-20! 

I_ 


V 3 “l_ 


1 

I 

i 

1 


i 

i 

i o 

i 



= ot 



where 0! and ot, can be found as follows 
o 1 

Eigenvalues of A are 

D 1 



= -1, X 2 = -2 ; therefore 

(X n ) k-1 = (_i) k_1 = a + a. X. = ot - ot 

v 1/ \ / o 1 1 o 1 



(X 2 ) k - 1 = (-2)*- 1 - o o + ^2 = a o -2^ 
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Substituting for k = 1,2, ... 



H x = 2 ( - 1) 1 1 - (- 2) 1 1 = 1 



2-1 2-1 

H 2 = 2 ( -1; 1 - (-2) Z 1 = 0 



2-1 3-1 

H 3 = 2 (-1) J - ( -2) J = -2 



which are exactly the same as the starting H^, k = 1,2, 
A second set of P and Q matrices is 





1 


0 




i 


0 


P 2 = 


0 


-h_ 


II 

CN 

a 


0 


1 



which give the realization 



x(k+l) = 



yOO = [l 



0 


-2 


x(k) + 


1 


1 


-3 




0 _ 










I — 1 


0 J 


l&M 





u(k) 



this realization also satisfies (17) . 

For the isomorphism of two realizations, W is calculated as 



w = (v 2 v 2 t ) v 2Vi t 
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By (13) 
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The two realizations have the following signal flow 



graphs : 




Y (z ) 



Y (z) 



Fig. 2 

»’ ' 

By Mason's Gain Formula both have the transfer function 

ti[z) (z+1) (z+2) 

By a long division 

-1 -2 -2 -4 -S 

H(z) = z + Oz - 2z + 6z - 14z ... 

the coefficients of z , k = 1,2, ... are again equal to 

Hj^, k = 1,2, ... 
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Then 



2. Example 


2 


Two— Inputs y 


Three-Outputs 








Continuous System 


(m = 2, 


p = 3) 


Let the 


V 


k = 1,2, 


t • • • 


be 








1 


0 




0 


1 




1 1 


H 1 = 


1 


2 


' H 2 " 


-2 


-4 


j ' H 3 * 


4 8 




2 




1 




-3 


-4 




8 14 




















-3 


-5^ 




7 


13 






H 4 " 


-8 


-16 


' H 5 " 


16 


32 








-18 


-34 




38 


74 






* 


1 


0 












HO. = 


1 


2 


P 1 


= 2 










2 


1_ 














1 


0 


0 . 


1 










1 


2 


-2 


-4 








H0 2 = 


2 


1 


-3 


-4 


p 2 


= 3 






0 


1 


1 


1 










-2 


-4 


4 


8 










-3 


-4 


8 


14_ 
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} •' 


... 1 


0 


0 


1 


1 


1 




1 


2 


-2 


-4 


4 


8 




2 


1 


-3 


-4 


8 


14 




0 


1 


1 


1 


-3 


-5 


H0 3 = 


-2 


-4 


4 


8 


-8 


-16 




-3 


-4 


8 


14 


-18 


-34 


* 


i' 


1 


-3 


-5 


7 


13 




4 


8 


-8 


-16 


16 


32 




1 

00 


14 


-18 


-34 


38 


74 


therefore 














n = 3 


and 


'q = 


2 










0 


1 


1 


1 








-2 


-4 


4 


8 








-3 


-4 


8 


14 






hot 2 = 


1 


1 


-3 


-5 








4 


8 


-8 


-16 








8 


14 


-18 


-34 







P. = 3 
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A set of P and Q matrices is 




1 0 



3 



-4 



Q 



1 



0 

0 



1 -2 

0 5 



3 

-7 



0 



0 -3 




which give the realization 

















— 






0 


1 


0 




l 


0 


• 

X 


= 


l 


1 


1 


X + 


0 


1 






-4 


-7 


-4 




0 


0 






1 


0 


0 








y 


= 


1 


2 


1 


X 










2 


1 


1_ 
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this has the impulse response matrix 



G ( t) 



3/2 - c + -| c 2t 



- 2 1 



-t -2 1 

2-3 € + 2e 



2C 



-2 1 



3/2 -2 c b + 5/2 c 2t 2-6 c + 5c 2t 



Expanding the exponentials in power series and expressing 
the result as a matrix polynomial the following expression 
can be obtained. 



G(t) 



1 


0 




0 


1 




1 1 


1 


2 


+ 


-2 


-4 


t + 


4 8 


2 


1 




-3 


-4 




i — 1 

“l 


-3 


-5 






7 


13 




-8 


-16 




+ 


16 


32 


fc4 /4I + * 


-18 


-34 






38 


74_ 





fc2 /2 



Each element of the coefficient matrices corresponds to the 
' s of (15) and the matrices themselves are exactly H^’s. 
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A second set of P and Q matrices is 




P 



2 



0 

-1 

2 

-2 




0 

0 

0 

1 

0 

0 



0 

0 

1 

-2 

2 

1 



0 0 0 

10 0 
0 0 0 

0 1-1 
2 10 
-1 -1 1 



Q 

2 




-3 

3 

-5 



3 

-2 

5 



3 



3 -3 



-4 

3 
-7 

4 



which give the realization 
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It can be shown as before that this realization also 

t 

satisfies (17) 



To show the isomorphism of two realizations W is 
calculated as 



V 1 “ C C 1 I A 1 TC 1 T ] • 





1 


1 i 


2 


0 


-2 


-3 


= 


0 


2 


1 


1 


-4 


-4 




_0 


1 


1 


0 


-2 


-3 




1 


0 


1 


0 


0 


0 


V 2 “ 


0 


1 


0 


1 


-2 


-1 




0 


1 


1 


0 


-2 


-3 








41 




8 


- 7 


‘ V 2 V 2 1 


1 -1 
) 


rH | t ^ 

II 


8 




29 


-16 








-7 




-16 


14 



V 2 V 1 



T 



3 

8 

16 



1 

15 

23 



1 

8 

15 
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w 



= (V 2 v 2 - 



-1 



V 2 V l' 



1 



0 



0 



0 

1 



1 0 
1 1 
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The W matrix also relates the two sets of P and Q 
matrices as 

















-i n 




W 


0 








w 


0 


P 2 = 


0 


I 


P 1 ' 


Q ; 


2 = Q 1 


0 


I 


this specific example 












1 


0 


0 


0 


0 


0 






0 


1 


0 


0 


0 


0 




P 2 = 


1 


1 


1 


0 


0 


0 


P 1 




0 


0 


0 


1 


0 


0 




*. 


0 

[ 


* 0 


0 


0 


1 


0 




v* 


r o > 


0 


0 


0 


0 


1 






1 


0 


0 


0 










0 


1 


0 


0 








Q 1 


-1 


-1 


1 


0 










_ 0 


0 


0 


1_ 









D. ANOTHER -ALGORITHM 

The algorithm given in [5] for discrete systems 
(equally applicable to continuous systems) with scalar 
measurement (single— output) shows a similarity with Ho ' s 
Algorithm. 
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For the algorithm, a free dynamical system (no forc- 
ing function) represented as in (18) is considered 
x(k+l) = A x(k) 

'v JJ ^ 

(18) 

y(k) = C x (k) 

^ u 



where 




a d 


is’ a nxn matrix 


C D 


is a lxn row vector 


X 


is a nxl state vector 


y 


is the scalar output 


and x(0) 


is known. Then HO, and HOT, matrices are construe 

i x, 



ted as in Ho ' s Algorithm for single-input^ single-output 
systems, from the measured output— sequence . An algebraic 
equivalent of the system is realized by letting 



a d : 
C D ■ 


= (HOT g ) (HO ) -1 

(19) 

= [l 0 o] 



Note that for single-input ^single-output systems of order 
n, q is equal to n, and p also is equal to n. Then this 
algorithm becomes a special case of Ho's Algorithm by let 



ting 

P = 
Q = 


I 

, (20) 
(HO ) ± 

q 



Because a nonsingular matrix can always be reduced to an 
identity matrix by elementary transformations on its 
columns (or rows) only. 



27 



,, For complete identification of the system, complete 
controllability , and observability was a necessary and suf- 
ficient condition for systems with forcing. Here, as 
is absent, controllability cannot be checked; instead a 
new condition called n-identification is imposed. 

A system is said to be n-identifiable if the matrix 

[x(o) | a d x(o) : a d 2 x(o) : ... ; A D n_1 x(o)] (21) 

has rank n. Similarity of (3) and (21) is obvious for 
single-input systems. 

I | 

Therefore it is possible to determine A and C if 

■ r * f < , . f JJ J-) 

the system is observable and n-identifiable. These two 

x. ; 

conditions are summed together and called-'l-identif iabili ty 
[5]. 



In fact another realization is possible by letting 



A" = (HO ) 1 (HOT ) 

D q q 

c" = [first row of HO ] 
D q 



( 22 ) 



This is also a special case of Ho ' s Algorithm; simply 



P = (HO ) 
Q = i 



(23) 



The following relation 





T 



( 24 ) 



I 
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holds for these two extreme cases. As HO and HOT are sym- 
metric matrices for single-input , single-output systems it 
can be written 



< a d )T = [< H ° T q > < h V" 1 ] T 
= [(no ,)' 1 ] 1 [ H 0 T q] T 

= (HO ) _1 (HOT ) = a" 
q 7 q 7 D 

Example 1 is applicable to these cases j 
letting 

x(°) = B d 



( 25 ) 



t simply 
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II. APPLICATION OF Ho 1 s ALGORITHM 



A. IMPLEMENTATION OF Ho ' s ALGORITHM 

Ho 1 s Algorithm perhaps is the one which requires min- 
imum computation among the existing identification schemes. 
Because of the simplicity of the computations it can be 
used in many applications. 

The algorithm is implemented with a digital computer 
program to see how it can be used in an off line fashion 
for identification. The results obtained can be extended 

for use in other applications. 

. . 5 

1. Digital Computer Programs 

The program consists of one main program and four 
subroutines. The main program reads the data, builds the 
HO and HOT matrices and controls the subroutines. For dig- 
ital simulation of systems, data can be generated within 
the main program by replacing the data reading loop with a 
generating function. The algorithm does not stop when two 
successive ranks of HO are equal, but uses all available 
data for additional checking. 

Subroutine SDWFD2 finds the successive derivatives 
evaluated at zero time for continuous systems; for dis- 
crete systems it is by-passed. FLAG controls this opera- 
tion. 



5 

See pages 56-61 
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Four Newton Forward Differences are used to evaluate each 
derivative. Subroutine RANK finds the rank of HO each 
time it is constructed. Subroutine PHOQ finds the elemen- 
tary transformation matrices P and Q which put HO into its 
normal form. Subroutine MULT makes the matrix multiplica- 
tions required by the algorithm. 

The program is written for single input-output systems, 
but can be altered for multiple input-output systems; in 
fact only the reading data and building HO and HOT parts 
need changes in the main program. No change is required 
in the subroutines. 

B. APPLICATION TO CONTINUOUS SYSTEMS 

For discrete systems the algorithm is straight forward, 
the only limiting factor is the accuracy of the measure- 
ments. For continuous systems however some difficulties 
arise. First, as a test input signal an impulse is needed. 
For linear systems the impulse response is the derivative 
of the response to a step input. Therefore a step input can 
be used as the test input. Higher order inputs can also be 
used if it is necessary as long as the right derivatives 
are obtained. The program is written to accept the step 
response of the system. 

A second problem comes from the numerical evaluation 
of successive derivatives from the measured samples. 
Numerical evaluation of a derivative is an approximation 
even with exact measurements, when higher derivatives are 
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needed accuracy diminishes greatly. A small error in the 
input data deteriorates the results, as the numerical meth- 
ods magnify these errors enormously. The sampling period 
is the biggest factor in this magnification as it appears 
in the denominator, and for each additional higher deriva- 
tive its power is increased by one. Magnification of error 
increases as sampling period decreases; but on the other 
hand for fast changing functions a small sampling period 
is needed. 

The examples presented in the following tables were 
simulated digitally for unit step-inputs to the systems. 

Two parameters, DX and DD, in addition to the sampling 
period, seemed to have great effect on the accuracy of the 
obtained results. A suitable choice of the three gave re- 
peating answers when the dimension of the HO matrix became 

— » r • 

larger than the order of the system. DX and DD are posi- 
tive small numbers; when finding the rank of HO, any 
number which has an absolute value less than DX is equated 
to zero by the program. When finding the P and Q matrices, 

which put HO to its normal form, any number which has an 

{ 

absolute value less than DD is made equal to zero by the 
program. The values of these parameters which gave a good 
result for this specific example are included in the 
tables. 

In example 3 a system with two real poles is consider- 
ed. The poles are not greatly separated and a good real- 
ization became possible because the step response is a 
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fairly smooth function, and in the evaluation of the deri- 
vatives not too much error is involved, in fact the first 
four derivatives are almost exact. In Example 4 a system 
with three real poles, which are separated by a significant 
amount, is considered. In part A, the sampling rate is not 
so high. Identification of the pole at the origin is 
almost exact, the pole at -1000 is identified fairly well 
but the far pole is in error by a large amount. In part B, 
the sampling rate is increased; while the pole at -1000 is 
identified more accurately and the far pole identification 
showed improvement, • identification of the pole at the ori- 
gin began to deteriorate. In Example 5, a highly damped 
complex pair of poles gave a very accurate identification 
with a high sampling rate; as there is no spread of poles 
a good choice of sampling rate results in an accurate real- 
ization. In Example 6, a very lightly damped complex pole 
pair is considered. Again identification is very accurate 
because of the same reasons in Example 5. 

In the above examples simulation of systems was carried 
out digitally in double-precision. In general a high samp- 

c. 

ling rate with a high-precision measurement is necessary 
which makes the scheme impractical to use with measurement 
of analog signals. A different approach is necessary which 
is the subject of the next section. 
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_ — 

EXAMPLE 3 


Transfer 
Func tion 
of The 
Sys tem 


300 

(s+5) (s+20) 


Sampling Per. 


0.005 


DX 


0.005 


DD 


0.000001 


Numerically 
Evalua ted 
Successive 
Derivatives 


0 . 3376xl0 -3 0.3000X10 3 -0.7500xl0 4 

0.1575xl0 6 -0.3187xl0 7 0.6393xl0 8 

-0.1278xl0 10 0.2448X10 11 0.2830xl0 13 


Realized 
Sys tem 
Matrix 


-0 . 2132xl0 -13 -0.1000xl0 3 

0.1000X10 1 -0.2500xl0 2 


Realized 
Dis tribution 
t Matrix ^ 


O.lOOOxlO 1 0.2541xl0“ 20 


Realized 

Observation 

Matrix 


0 . 3376xl0 -3 0.3000X10 3 


Eigenvalues 
of System 
Matrix 


-0.4999 x 10 1 + j 0.0 
2 

-0.1999 x 10 + j 0.0 
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EXAMPLE 4a 


j 

Transfer 
Function 
of The 
System 

t 


io 12 

s (s+1000) (s+1000000) 


i 

i 

Sampling Per. 


0.001 


DX 


1.0 


DD 


1.0/| largest element of HO | x 2.0 


j 

Numerically 

Evaluated 

Successive 

Derivatives 


0.430xl0 2 0.9148xl0 6 0.8746xl0 9 

0.8361xl0 12 -0.7993xl0 15 0.7641xl0 18 

21 24 27 

-0.7304x10 0.6981x10 -0.6670x10 


j 

Realized 

System 

Matrix 


-0 . 5684xl0 -13 -0 . 1309xl0 -9 0 . 2131xl0 -1 

O.lOOOxlO 1 0.3725xl0 _8 -0.1991xl0 7 

0 . 1065xl0 -13 O.lOOOxlO 1 -0.3039xl0 4 


Realized 
Dis tribution 
Matrix T 


O.lOOOxlO 1 0 . 3552xl0 -14 -0 . 1387xl0 _16 


Realized 

Observation 

Matrix 


0.430xl0 2 0.9148xl0 6 -0.8746xl0 9 


Eigenvalues 
of System 
Matrix 


-0.2083 x 10 4 + j 0.0 
-0.956 x 10 3 + j 0.0 
0.1070 x 10 _7 + j 0.0 
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EXAMPLE 4b 


Transfer 
Func tion 
of The 
System 


io 12 

s (s+1000) (s+1000000) 


Sampling Per. 


0.0001 


DX 


O 

• 

i — 1 


DD 


1.0/| Largest element of HO | x2.0 


Numerically 

Evaluated 

Successive 

Derivatives 


-0.9625 O.lOOOxlO 7 -0.9918xl0 9 

0.8123xl0 12 0.2927xl0 16 -0.8084xl0 20 

0.1704xl0 25 -0.3552xl0 29 0.7400xl0 33 


Realized 

System 

Matrix 


0 . 863 3x10 -12 -0 . 5 544x10 -8 -0.5187xl0 3 

O.lOOOxlO 1 -0 . 19 89x10 _11 -0.2083xl0 8 

-19 1 5 

0.6098x10 0.1000x10 -0.2183x10 


Realized 
Distribution 
Matrix T 


O.lOOOxlO 1 -0 . 596 3x10 ~ 18 -0 . 2646xl0 _22 


Realized 

Observation 

Matrix 


-0.9625 O.lOOOxlO 7 -0.9918xl0 9 


Eigenvalues 
of System 
Matrix 


-0.2083 x 10 5 + j 0.0 
-0.9999 x 10 3 + j 0.0 
-0.2490 x 10~ 4 + j 0.0 



36 





EXAMPLE 5 


Transfer 
Function 
of The 
System 


1000000 

(s+25000 + j 25000) (s+25000 - j25000) 


Sampling Per. 


0.000001 


DX 


0.000005 


DD 


1.0/| Largest element of HO | x 2.0 


Numerically 

Evaluated 

Successive 

Derivatives 


0.6764X10 1 0 . 8164xl0 6 -0.4927X10 11 

0.1443xl0 6 -0.1056xl0 20 -0.1306xl0 25 

3 37 44 

0.3408x10 -0.2533x10 0.2595x10 


Realized 
Sys tern 
Matrix 


0 . 9094xl0 -11 -0.1250xl0 10 

O.lOOOxlO 1 -0 . 5000x10 5 


Realized 
Distribution 
Matrix T 


1 -20 
0.1000x10 0.1058x10 


Realized 

Observation 

Matrix 


0.6764X10 1 0.8164xl0 6 


Eigenvalues 
of System 
Matrix 


-0.2500 x 10 5 + j 0.2500 x 10 5 
-0.2500 x 10 5 - j 0.2500 x 10 5 
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EXAMPLE 6 


Transfer 
Function 
of The 
Sys tem 


1001 

(s + 0.01 + j 500) (s + 0.01 - j 500) 


Sampling Per. 


0.0001 


DX 


0.005 


DD 


0.000001 


Numerically 

Evaluated 

Successive 

Derivatives 


-0 . 2084xl0~ 6 -O.lOOlxlO 4 -0.1986xl0 2 

-0 . 2502xl0 9 0.1128xl0 8 0.6234xl0 14 

0.2890xl0 17 -0.3539xl0 22 0.4010xl0 27 


Realized 

System 

Matrix 


0 . 8198xl0 -13 -0.2500xl0 6 

0.1000X10 1 -0 . 1984xl0 _1 


Realized 
Distribution 
Matrix ^ 


O.lOOOxlO 1 -0 . 2082xl0 -9 


Realized 

Observation 

Matrix 


-0 . 2084xl0~ 6 O.lOOlxlO 4 


Eigenvalues 
of System 
Matrix 


-0.9922 x 10“ 2 + j 0.4999 x 10 3 
-0.9922 x 10" 2 - j 0.4999 x 10 3 
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C. DISCRETE REALIZATION OF CONTINUOUS SYSTEMS 

As is seen from the preceding examples, to find a set 
of differential equations for continuous systems is some- 
what a trial and error procedure, . as one must come up with 
a suitable triple of DX, DD and sampling period, if no 
apriori knowledge is available for the system. More impor- 
tantly, very accurate measurements are needed which is im- 
possible in an actual experiment with hardware. 

The main error source in the above procedure is the 
evaluation of derivatives numerically. If the measured 
data can be used directly, without evaluating derivatives, 
inclusion of this error to the algorithm can be. eliminated. 

This approach leads to a system which is operating con 
tinuously but observed discretely as in Fig. 3 




y(kT) 



u(t) 


Continuous 








System 





y (t) 



Fig. 3 

where T is the sampling period and u(t) is a step. If the 
system has the dynamics as in (26) 

x = Ax + Bu 

(26) 

y = Cx 
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Then 



x(t) = X {(sI-A) 1 x(t Q ) + (sI-A) 1 Bu ( s ) | (27) 



As u(t) is a step (amplitude V) 



x(t) -X { ( s I - A ) 1 x(t o ) + ^ Sl s A) -v} 


(28) 


Defining 

0 (t-t Q .) =oL- • \ 


[(sI-A) -1 } 




r(t-t o ) ^JL" 1 \ 




(29) 


Substituting in 






r x(t) = 0(t-t )x(t )+ F(t-t ) V 
^ o ^ o o 


(30) 


letting 






t = kT 
o 

t = t + T 
o 




(31) 


x[ (k+1) t] = 0 (T) x (kT) + F(T)V 

~ r\_/ 


(32) 


and 






y (kT) = Cx (kT) 




(33) 



0 (T) and F(T) are constant matrices as the sampling 
period is constant, V can be replaced by v(kT) in (32) let- 
ting 

v(kT) =1 , k = 0,1,2, ... (34) 

as the test input is a unit step. 
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Dynamics described by (32) and (33) can be represented 



as a discrete system of Fig. 4. 



v(kT) 



Discrete 



y(kT) 



■> 



Sys tern 



Fig. 4 



If Fig. 3 is driven by a unit step and Fig. 4 with a 
unit impulse train of period T, they both will give the 
same y(kT) . Therefore samples taken as the result of ex- 
periment on Fig. 3 can be used to identify the system in 
Fig. 4 as if it is driven by a unit impulse train and the 
output is measured. 

Dynamics described in (32) and (33) is an exact re- 
presentation of the continuous system of Fig. 3 at sampling 
instants and constitutes a discrete realization of the con- 
tinuous system. Two points of importance must be mentioned 
here. One is the loss of controllability and observability 
due to the sampling of continuous systems [2], as a result 
a realization which is not a faithful representation of 
the system will be obtained. If a linear— stationary— con- 
tinuous system is controllable and observable, it will re- 
main so when sampling is introduced if and only if [6] 




k2 n 



( 35 ) 



T 
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where 



s , m = 1,2, ...n eigenvalues of system 

m 

i / j = 1/2, n 

k = positive integer 
T = sampling period 

If no apriori knowledge is available about the system, at 
least two realizations must be obtained with different 
sampling rates to avoid misrepresentation. 

The second point comes from the usage of a unit step 
as the test input. As a result the discrete equivalent in 
Fig. 4 is driven with. a unit impulse train. If the ori- 
ginal system has a transfer function 



H(z) = = a ± z 1 + a 2 z 2 + a 3 z 3 + ••• (36) 

in the z-domain, the realization will be for 



h' (z) 



zN(z) 

(z-l)D(z) 



a i Z 



-1 



. » -2 
+ a 2 z 



i -3 
+ a^z 



+ 



(37) 



a^ , i = 1,2, ... are measured sample values to a unit 

J 

step input. This will give an additional eigenvalue in- 
creasing the order of the system by one. This can be 
avoided by a simple manipulation of the data. 



H(z) = h' (z) (1-z ■*■) = a|z + (a 2 ~a|)z 2 + ••• (38) 

Therefore 

I 
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This manipulation of the data has an additional benefit, if 
there is an unknown biasing in the measurement it is remov- 
ed automatically except from the first sample. The program 
takes care of this conversion of the data. 

In the following examples systems are simulated digit- 
ally. To see the sensitivity of the algorithm to numerical 
precision, sample values were rounded-off to four signific- 
ant digits, and a realization without manipulation of the 
data is included in some examples to make comparisons pos- 
sible. As will be seen the choices of DX and DD are not 
critical as in the continuous representation of systems. 

Example 7a shows the identification of the additional 
pole at the origin because of a unit step test signal. In 
7b this pole is removed with the manipulation of the data, 
and six significant digits are used for sample values, 
realization is very accurate. In 7c sample values are 
rounded to four significant digits, which results in a de- 
terioration in the answer. In Example 8a and 8b seven 
significant digits are used. In 8a an additional pole is 
present because of step input, in 8b it is removed; with 
the same significant digits 8b gave a more accurate answer 
as the system order decreased by one, less computation is 
required and error is decreased. In 8c four significant 
digits are used with some deviation of the answer as a re- 
sult. In Example 9a and 9b a complex pole pair and a real 



43 



pole are considered, with eight and four significant digits 
of sample values respectively. Highly accurate results are 
obtained in both cases. Example 10a shows the loss of con- 
trollability and observability (as the sampling is intro- 
duced in the measurement, only observability is lost) due 
to the sampling. In 10b with a different sampling rate the 
system is identified completely. 
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EXAMPLE 7 a 


Transfer 
Function 
of The 
Sys tem 


160 

(s+1) (s+2) 


Sampling Per. 


0.1 


DX 


0.001 


DD 


0.00001 


The Type 


Not manipulated to take the step out 


of The Data 


Double-Precision 


Realized 
Sys tem 
Matrix 


0 . 3552xl0 -14 0.0 0.7408 

0.1000X10 1 0 . 3552xl0 -14 -0.2464X10 1 

0 . 6661xl0~ 15 0.1000X10 1 0.2723X10 1 


Realized 
Distribution 
Matrix T 


0.1000X10 1 0.3552xl0" 14 -0 . 2220X10 -35 


Realized 


I i 


Observation 

Matrix 


0.7244 0.2628x10 0.5374x10 


Eigenvalues 


0.8187 + j 0.1230 x 10 -14 


of System 


0.9048 - j 0.7288 x 10 _15 


Matrix 






0.1000 x 10 1 - j 0.2778 x 10 _15 


Corresponding 


-2.0 


s -Domain 
Eigenvalues 


-1.0 




0.0 
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EXAMPLE 7b 


Transfer 
Function 
of The 
System 


160 

(s+1) (s+2) 


Sampling Per. 


0.1 


DX 


0.001 


DD 


0.00001 


The Type 
of The Data 


Manipulated to take the step out 
Rounded to six significant digits 


Realized 

System 

Matrix 


0.2220 x 10~ 15 -0.7408 

0.1000 x 10 1 0.1723 x 10 1 


Realized 
Distribution 
Matrix T 


-0.-1000 x 10 1 -0.5551 x 10 -16 


Realized 

Observation 

Matrix 


0.7244 0.1904 x 10 1 


Eigenvalues 
of System 
Matrix 


0.8187 + j 0.0 
0.9048 + j 0.0 


Corresponding 
s -Domain 
Eigenvalues 


-2.0 

-1.0 
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EXAMPLE 7c 


Transfer 
Func tion 
of The 
System 


160 

(s+l) (s+2) 


Sampling Per. 


0.1 


DX 


0.001 


DD 


0.00001 


The Type 
of The Data 


Manipulated to take the step out 
Rounded to four significant digits 


Realized 
Sys tem 
Matrix 


0.6661 x 10 -15 -0.7386 

0.1000 x 10 1 0.1722 x 10 1 


Realized 
Dis tribution 
Matrix T 


0.1000 x 10 1 -0.6938 x 10 -16 


Realized 

Observation 

Matrix 


0.7245 0.1904 x 10 1 


Eigenvalues 
of System 
Matrix 


0.8069 + j 0.0 
0.9154 + j 0.0 


Corresponding 
s -Domain 
Eigenvalues 


-2.145 

-0.885 
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EXAMPLE 8a 


Transfer 
Function 
of The 
System 


4000 

(s+1) (s+5) (s+10) 


Sampling Per. 


i — 1 

• 

o 


DX 


0.001 


DD 


0.00001 


The Type 
of The Data 


Not manipulated to take the step out 
Rounded to seven significant digits 


Realized 
Sys tern 
Matrix 

' ; 


0.426xl0" 13 -0.5684xl0 -13 -0 . 5684xl0 13 -0 . 2030 
0.1000X10 1 -0.5684xl0 13 -0.5684xl0 13 0.1310X10 1 
0.0 O.lOOOxlO 1 0.5684x10 13 -0.2987x10 1 

-0.3907x.10 -13 0.0 O.lOOOxlO 1 0.2880X10 1 


Realized 
Distribution 
Matrix T 


-) 

0.1000X10 1 0.0 0.0 0.0 


Realized 

Observation 

Matrix 


0.4537 0.2542X10 1 0.6169X10 1 0.1077xl0 2 

, f 


Eigenvalues 
of System 
Matrix 


0.9987 - j 0.4198 x 10 -14 
0.9073 - j 0.2743 x 10 -13 
0.3715 - j 0.2742 x 10 _13 
0.6028 + j 0.5936 x 10 _13 


Correspond- 

ing 

s-Domain 

Eigenvalues 


0.0 

-0.972 

-9.9 

-5.06 
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EXAMPLE 8b 


Transfer 
Function 
of The 
Sys tern 


4000 

(s+1) (s+5) (s+10) 


Sampling Per. 


i — 1 
• 

o 


DX 


0.001 


DD 


0.00001 


The Type 
of The Data 


Manipulated to take the step out 
Rounded to seven significant digits 


Realized 

System 

Matrix 


0.2220xl0 -15 -0.1554xl0" 14 0.2018 

0.1000X10 1 0.3552xl0~ 14 -0.1104X10 1 

0.2220xl0 -15 O.lOOOxlO 1 0.1879X10 1 


Realized 
Distribution 
Matrix T 


0.1000X10 1 0.1554xl0 -14 -0 .4440xl0~ 15 


Realized 

Observation 

Matrix 


0.4537 0 . 2088x10 1 0.3627X10 1 


Eigenvalues 
of System 
Matrix 


0.9047 + j 0.3596 x 10~ 12 
0.3676 - j 0.3226 x 10~ 12 
0.6067 - j 0.3617 x 10 -13 


Corresponding 
s -Domain 
Eigenvalues 


-1.001 
-lb. 00 3 
-4.999 
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EXAMPLE 8c 


Transfer 
Function 
of The 
Sys tem 


4000 

(s+1) (s+5) (s+10) 


Sampling Per. 


0.1 


DX 


0.01 


DD 


0.00001 


The Type 
of The Data 


Manipulated to take the step out 
Rounded to four significant digits 


Realized 
Sys tem 
Matrix 


0 . 133 2xl0~ 14 0 . 1776xl0 -14 0.1896 

O.lOOOxlO 1 0.0 -0.1080X10 1 

r 

0 . 1554xl0 -14 O.lOOOxlO 1 0.1866X10 1 


Realized 

Dis tribution 
Matrix T 


O.lOOOxlO 1 -0 . UlOxlO -14 -0 . 2220xl0 -15 


Realized 

Observation 

Matrix 


0.4537 0.2088X10 1 0.3628X10 1 


- Eigenvalues 
of System 
Matrix 


0.8977 - j 0.1602 x 10 _14 
0.3316 - j 0.1348 x 10 _14 
0.6368 + j 0.3664 x 10 _14 


Corresponding 
s -Domain 
Eigenvalues 


-1.08 

-11.0 

-4.51 
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EXAMPLE 9a 


Transfer 
Function 
of The 
Sys tem 


80000 

(s+40) (s 2 +4s+400) 


Sampling Per. 


0.1 


DX 


0.01 


DD 


0.00001 


The Type 
of The Data 


Manipulated to take the step out 
Rounded to eight significant digits 


Realized 
Sys tem 
Matrix 


0.1942xl0~ 15 -0. 4440x10 _15 0.1227xl0 -1 

O.lOOOxlO 1 -0 ,4440xl0 -15 -0.6581 

0.0 O.lOOOxlO 1 -0.6481 


Realized 
Dis tribution 
Matrix T 


O.lOOOxlO 1 0 . 1332xl0 -14 0.1332xl0 -14 


Realized 

Observation 

Matrix 


0.4355X10 1 0 . 3682x10 1 -0.4645X10 1 


Eigenvalues 
of System 
Matrix 


-0.3332 + j 0.7478 
-0.3332 - j 0.7478 
0.1831 x 10 _1 + j 0.3330 x 10 -15 


Corresponding 
s -Domain 
Eigenvalues 


-1.998 + j 19.9 
-1.998 - j 19.9 
-40.0 
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EXAMPLE 9b 


Transfer 
Function 
of The 
System 


80000 

(s+40) (s^+4s+400) 


Sampling Per. 


0.1 


DX 


0.01 


DD 


0.00001 


The Type 
of The Data 


Manipulated to take the step out 
Rounded to four significant digits 


Realized 
Sys tem 
Matrix 


-0 . 888 1x10 -1 5 -0 .4440x10 -1 5 0.1263xl0 _1 

O.lOOOxlO 1 0.0 -0.6576 

-O.lllOxlO -14 O.lOOOxlO 1 -0.6474 


Realized ■, 
Distribution 
Matrix T 


O.lOOOxlO 1 0 . 6661xl0 -15 0 . 2220xl0 -15 


Realized 

Observation 

Matrix 


0.4356X10 1 0.3683X10 1 -0.4646X10 1 


Eigenvalues 
of System 
Matrix 


-0.3331 + j 0.7478 
-0.3331 - j 0.7478 
0.1885 x 10 -1 + j 0.5551 x 10 -16 


Corresponding 
s -Domain 
Eigenvalues 


-1.998 + j 19.9 
-1.998 - j 19.9 
-39.8 
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EXAMPLE 10a 


Transfer 
Function 
of The 
System 


80000 

(s+40) (s+2+jlOX) (s+2-jlOX) 


Sampling Per. 


0.1 


DX 


0.00001 


DD 


0.00001 


The Type 
of The Data 


Manipulated to take the step out 
Double-Precision 


Realized 
Sys tern 
Matrix 


-0.2220 x 10 -15 0.1499 x 10 -1 

0.1000 x 10 1 -0.8004 


Realized 
Distribution 
Matrix T 


0.1000 x 10 1 -0.2220 x 10 _15 


Realized 

Observation 

Matrix 


0.2981 x 10 1 -0.1765 x 10 1 


Eigenvalues 
1 of System 
Matrix 


-0.8187 + j 0.0 
0.1831 x 10 -1 + j 0.0 


Corresponding 
s -Domain 
Eigenvalues 


-2.0 + j 10 X 
-40.0 
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EXAMPLE 10b 


Transfer 
Function 
of The 
Sys tem 


80000 

(s+40) (s+2+j 10X) (s+2- j 10X) 


Sampling Per. 


0.05 


DX 


0.00001 


DD 


0.00001 


The Type 
of The Data 


Manipulated to take the step out 
Double-Precision 


Realized 

System 

Matrix 


-0 . 8881xl0~ 15 -0.2775xl0 -16 0.1108 

0.1000X10 1 0.1387xl0 -16 -0.8187 

. -0.4440xl0" 15 O.lOOOxlO 1 0.1353 


Realized 
Distribution 
Matrix T 


O.lOOOxlO 1 0.1387xl0~ 16 0.5551xl0 15 


Realized 
Observation , 
Matrix 


0.8901 0.2091X10 1 -0.1333 


Eigenvalues 
of System 
Matrix 


0.8471 x 10 -7 + j 0.9048 

0.8471 x 10“ 7 - j 0.9048 

-14 

0.1353 - j 0.9048 x 10 


Corresponding 
s -Domain 
Eigenvalues 


—2.0 + j 10X 
-2.0 - j 10X 
-40.0 
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III. CONCLUSION 



Ho's Algorithm provides a valuable tool for complete 
identification of linear— stationary systems which are com- 
pletely controllable and observable. The computations re- 
quired are simple and can be programmed easily, in fact for 
most of the computations, available library subroutines of 
a computer facility can be used. 

For discrete systems accuracy of the measurement is 
the only limiting factor. For continuous systems a dis- 
crete realization seems to be more feasible as the evalua- 
tion of derivatives introduces error. The tested examples 
showed that a f ive~digits-precision data would give fairly 
good results for this approach. 
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onn ry n n;-> n n n no ~>on nnonn ^>r>n n 



DIGITAL rnwpMTFP RR0GRAMR 



m/\ t n r q nr, PA M 

pt*,n<; T|-ip p^T* PNJF n M FAOH CAPO w^TH 
F n P M tT(P17, 10) . FPP DIRCLFTF RE p Pr-S FNT IT ION OF 
TM r CYS T FM fi T L c AST ?*M SAM p LFS MUST R F 
Rt)P D L I FP p I pc T BFTMG RAM p L c D ONE FERIOC A FT FP 
INITIAL time. FOP continuous RFPRFS c NT AT ION 
PF THF cyjpM (EVALUATING O c R TV AT !V p S ) AT LF/ST 
° * M+ 1 SAMFIFS “UST B p SUPPLIED FIPST BEING 
CAMPLFD at INITIAL ttmf, WHF R c w AT ! FAST 
TMC c STI Mf. T C D UPPpR LIMTT HP SYSTEM HI M P N C . ION 
FLU C PNF. 

cppGRAM J~ FOP c I NO | F IMPIIT SINGLE OUTPUT 
SYSTcwf, pijt f,*N RT ALTERED F^d mijijIPLF 

T NptJT-nUTPlIT <CYST c MC. 

5 PP ROUTT NF S PO NOT Nfco ANY CHANGr. 



I M p L I P I T PEA L*p ( A-H, p -7 ) 

CO'Mdlc X*1 A A X , VA L 

nr MPMC TON MT AO) ,H p (2 p ,?n) ,HO p (2 0 ,20 ) , P ( 2C , 20 ) , 0 ( 20 , 20 ) 
fi »* A ( 20,20 , RBI? 0,20) ,C.o (20 ,2 0) , A A1 (20 ,20 ), HOI (20, 20 ), 
FA( l* c ) , A X ( 20, 20 , VA L ( 20 ) 

MM I p PIM^MCION OF I N°UT VFCTOR 
MM= l 

mc'ic DIMFNSIPN OF OUTPUT VECTOR 

MP= 1 

FLAG CrNTPPl c MQPP OF ALGORITHM, p LAG=l.O GIVES A 
p FT n p niFFFRfNCE FOUATT 0 N C , FLA 0 ,= 0.0 GIVES A $E T 
pr ni c Fcpc N xiAL EQUATIONS. 

P L A 0 = 1 . 0 

r m if TH C FSTIMATFO i i° f c R LIMIT OF T H c SYSTEM hrqfr 



c 


p L U * 


ON? AT 


[FAST 




M= c 








| =M*2 
L_n = 4*L + 1 






r 


W T * 


ffi mpli MG PERIOD 




W = 0. 1 
WR I tf < a , 


20) V, 






TF( FLAG . 


GT.C.O) 


on to ci 


c 


f f fir 


r A M p LE 


Vf l UFS 




on o I = l 


tin 





pEAP/( c , 111) A(I) 
a CPMT T NUF 

r F I Mp TUCCFSSIVE DFPIVATTVFP 

CALL SO WFP2 I H , W , L , A ) 

GO Tr <=3 

r r f a o sample value* 

C 1 DO 12 1 = 1 ,L 

ROA0( c ,lll) A ( T ) 

12 OONTINUF 

C FIND P ULSF RFSPPMRF S AMPL P VALUE* 

H( 1 )=A ( 1 ) 
n 0 02 I. L = 2 , L 
c 2 H(|.L) =A (LL)-A(LL-l > 

S3 VP I TF ( f , 400 ) 

WRTTP( Af AC) (Hf J) ,.J = 1 ,L) 

DO 1^ K=1,M 

C RUT l T HC AND HOT MAtrjcf* 

DO 14 1=1, K 

on 1 / .1 = 1 , k 
N 1 . = I + . I - 1 

HO T ( I , J)=H(NL+l) 
u 0 1 ( I , J ) = u ( NL) 

14 H0( I , J)=H( NL) 

K 1=K*MM 
D X = 0 . C C 1 
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dd=c. ococi 

C PIMP THP RANK of WO 

CALL RANK(K,K1 ,H0l ,CX,IC) 

WP I TF ( f , C C ) IC,K 
C F I MO P AND 0 MATRICES 

CALF DHrO«K f Kl f Hr t DD,PfO,IC) 

C MULTIPLY d f, mo Hnr 

CALI mui t, k,p,HC t ,K,K1 ,AA1 ) 

WP I tc ( f- ,A0) 

C multiply THF PF SUL t WITH q 

CALL MUL T < Kl.AA 1 ,0,IC ,IC ,AA) 
r WP I TF * YS TF M MATRIX 

no 4 I A = 1 , I C 

4 WR ! Tp ( * 1 1 Q ) ! A A ( I A , J 4 ) , JA = 1,IC) 

W R I TP ( 6, C C3) 

WP I t p ( 6,70 

C MULTIPLY F AND HP 

CALL MULTI K,P,Mr,IC ,MM t AR) 

C WP I TF DISTRIBUTION MATRIX 

nn 2 I A = 1 , I C 

2 WRI T F(Oin (BB(IA,JA) ,JA = l t MM) 

WPTTC( 6,500) 

WR ITC I 6, ?C) 

C MULTIPLY HO AND 0 

CALL mul t ( K1 ,H0 , 0 , MP , I C , CC ) 

C. WP I TF OBSERVATION MATRIX 

DO 3 !A=1,MP 

3 WRITE!*, 1C) (CC IT A , JA ) , JA=1 , IC ) 

FIND THE EIGENVALUES OF SYS T EM MATRIX. OALMAT IS A 
L T P R A P Y S UB R OUT I N c (N.D.G.S. COMPUTER FACILITY) 
WHICH FINDS THE COMPLEX EIGENVALUES OF 4 COMPLEX 
MA TRT X 

DO c T A = 1 , I C 
DO 5 JA = 1 , T C 

5 A X I T A , JA ) =AA 1 1 A ,JA) 

CALL DALMATI AX,VAL ,IC ,20,NCAL) 

WP ! T ( 6 , °0 ) IVAL(IA) ,IA=1 , IC) 

1 K CONTI MUF 
10 FnRMATi /,7D18.10) 

2C F0Rmat( ///T 2C, ' SAMPLING PE R I 0P= • , 01 6 .3 , IX , ’ S ECONDS • f / ) 
4C FORMAT! 7018. 10) 

K 0 RQRM.A T! //T4, 'RANK H0= • , I 3 , 5 X , ' DI MENS ION H0=',I3,///) 

6 C F0RMAT(TA t *4 MATRIX’,//) 
t 0 FORMA T(T4,*R MATRIX’,//) 

BO FflR MA T ( T4 , • C MATRIX’,//) 
cc FORMAT! //,10X, ’EIGENVALUE*’ ,2018.10) 

111 FORMA T<riT. IQ) 

400 form; T( //T20, »H VALUES’,/) 

*CC PORMA T( /////) 

RETURN 
PND 



SUBROUTINE SO WF 02 
PUR PfiPF : 

SUBROUTINE S0WF02 FINOS THE SUCCESSIVE 
DFRIVATIVFS OF A FUNCTION EVALUATED AT 7FR0 
VALUE OF THE I NOEPENOEN T VARIABLE U^ING 
NEWTON FORWARD DIFFERENCES. 

PESCPIPTTCN OF PARAMETERS: 

A — ORDINATE VALUES OF Fll NC T I ON 
SAMPLED WITH period W 

in — highest derivative t 0 be found 

W — SAMPLING PERIOD 

H — CONTAINS SUCCESSIVE DERIVATIVE* 

® , 1 0 , W A.RF I N C UTR ,H IS OUTPUT 



NOTE : 

FPR I DTH DERIVATIVE 4*ID+1 SAMPLES AR F NFEDED 
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subroutine sdwfd2(h,w,io,a) 

IMPLICIT REAL*°(A-H,°-2) 

DIMENSION HI 4()) ,A(16E),B(lA t >),Cll65),D(lA5),E(lAE) 

C IMTI.5LI7F THE INOFCTS 

K=4*yn 
|_ = 4*T n- l 
M = A*T f)_2 
M=4«I n -3 

on ? j=i,in 

c FIMO TH C F I P S T DIFFERENCES 

DC 2 I=1,K 

2 R( ! ) = A( I + D-M I ) 

c FIND TH C SET OND 01 F F p R E NIC F S 

on 3 I = 1 1 L 

3 Cl I )=B( I+1)-B(I > 

C FIND THE THIRD DI FFFR r NCES 

no 4 i =i , v 

4 0 ( I )=C( I+D-CIT ) 

C FIND THF FOURTH DIFFERENCE 

DO 5 I = 1 , M 

5 E ( T >=0( I + 1 ) -0 1 l ) 

CALCULATE THF Jth DFRIVATIV C AND PII T IT to CAMPLE 
AP PAY TO F I NO THE NEXT "PRIVATIVE 
no 6 I = 1,N 

a Am = i p r n-o.5*c m + ( i . 0 / 3.0 *m 1 )-o .25 *ei i > > /w 

PUT th 1 : DERIVATIVE FVAUIATEO AT 
7 FR C TO DERIVATIVE ARRAY 
HI J ) = A ( 1 ) 

K=K-4 
L = L-4 
M = M-A 
7 N=N-4 
RETURN 
FND 



SUBRCU T TNR RANK 



PURPOSE : 

SUBROUTINE RANK FINDS THE RANK OF REAL MATRICES 
WITH DIMENSIONS UPTO (20 ,20) 

OF SC R I RT I CN OF PAPAMFTFRS: 

A — INRUT MATRIX DIMENSIONED RFAL*P AI20,20) 

M — ACTUAL POW niMFNSIQN n<= A 

N — ACTUAL COLUMN DIMENSION OF A 

jr — A T t H p end ic is THE RANK OF A 

DD — A SUITABLE S M A L L POSITIVE NUMBER, DUR ING 

OPERATION ANY NIJMRER WHICH HAS AN ABSOLUTE 

V*LUF LESS THAN DO IS MADE EQUAl TO 7FR0 

A,m # n,DD APE IMPUTS IC IS OUTPUT 



NOTE : 

THF A MATRIX IS DESTROYED 



SUBROUTINE RANK ( M , N , A , DD , I C ) 

IMPLICI T REAL*fi(A-H,R-7) 

DIMENSION A I 20 , 20) 

IC INITIALIZED TQ 7FR0 AT THE START 
IC = 0 

DO 14 I P = 1 ,M 

FACH TIME BIG IS INITIALISED AS 
ftrct E L c ME NT OF TRTH ROW 
B I G=DA R SI A I I R , 1 ) ) 

I = I P 
J=1 

SFAPCH THE MATRIX FOR LARGEST ABSOLUTE VALUE 
ELEMFNT BELOW THE IRTH ROW (INCLUDING) 

DO S L = 1 , N 

DO 5 K=TR,M 
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IF( big.gf .oabs( c (K,L) ) ) go to «■ 

C BIG I S TH r LARGEST ABSQLUTF VALM C 

r rj cyrMT nuRING RF*PCH 

B I G = OA P S( A ( K , L ) ) 

C MfMPPf7F THF DpW 4 NP COLUMN OF PIG 

I=K 
)=t 

3 CDMTTNt.r 

C TP PIG IS IN T P TH ROW GO TO Q , OTHERWISE INTERCHANGE 

C THE pr, W OF RIG A NO THE IRTH ROW 

IP( T .FO. I® ) GO TP 1 P 

00 7 L = 1 » N 

TFMCV =A ( T ,u 

Af I ? t ) = A f I R » L ) 

7 A( !R , l )=Tr«oi 
R T F ( P IG . I T .nn| GO TO 1^ 

C IP BIG IS NOT 7ER0 !NCR C ASE RANK ONF 

ir = ir + 1 

r P F PUC p ALL P L r MR NT S ABQVP AMP BFLOW OF BIG T 0 ZERO 

nn 12 K=1 t w 

TF(K.Pn.iR) GO Tn 12 

1 F ( A ( K , J) . = 0.0.0) GO TO 12 
G = M K , J) / » ( IP, J) 

op 11 L = 1 , N 
IPfL.FQ.J) GO TO 13 
A(K,L)=A(K,L> — g *A ( I R , L ) 

IF(DAB C (A(K,L) ) .LT.DD) A(K,L)=0.0 

13 A ( K , J ) = C . 0 

11 CONTI M IJF 

12 CONTTNI'F 

14 CONTI NIJF 

p p T(JP N 

C N0 



SUB P OUT T MF PHOO 
ClioprRF : 

SUBROUTINE RHHQ FINOS THE ELEMENTARY 
transformation MATRIC C S o ano 0 WHICH 
PUT THE A MATRIX IN x O ITS NORMAL FORM. 

OESCPIPTTCN OF PA o A MPTF RS : 

A — THP MATRIX ITS NORMAL FORM WILL BF FOUND 

V ACTUAL ROW DTMPNSTON CP a 

N — ACTUAL COLUMN Dl mf NS ! CN OF A 
I C — RANK OF A 

p — F L F mp NT.A R Y orw TRANSFORMATION MATRIX 
0 — f l.F MF NTA R Y COLUMN TRANSFORMATION MATRIX 
on — f SUITARLF SMALL POSITIVF NU MPER , OUR I NG 
OPERATION AMY NUMBER WHICH HAS AN ABSOLUTE 
VALUE LESS THAN DO IS MADE EQUAL TO 7ER0 
M , N , A ,DD , IC APE INPUTS, P AND 0 A R F OUTOIJTS 

mote : 

A ,R AND 0 MATdjcES Apj: ORST ARRANGED IN THE X 
ARRAY SUCH TH A T O I s on T 0P AND P ON LEFT OF A. 
P ANn o ARE IDENTITY MAT R IC FS AT THE BEGINING. 



SUBROUTINE PH00(M,N,A ,DD,P,Q,IC) 

IMPLICIT REAL*9( A-H,P-7) 

DIMENSION A (20, 20 ,P (20, 20 ,0(20,20) ,X (40,40) 
C FIND THE DIMENSION OF X 

K 1=M+N 

C INITIALIZE X TO 7FR0 

HP 1 1 = 1 , Ki 

DO 1 J = 1 , K 1 
1 X( I , J ) =0 . 0 

C INITIALIZE P PART OF X TQ IDENTITY 

K 2 = N + 1 

Do 2 !=K2 ,K1 
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DP 2 J«1,M 

IF( J.NF.f T-N ) » GO TP 2 
X( I t J ) =1 . C 

2 CONTINUE 

r I N I T 1 4 1 1 7F 0 PART OF X TO IDENTITY 

K3=M+l 
^P 3 1=1, N 
0^ 3 J=K3,K1 
I F ( I , N r , ( J-M) ) r,0 TO 3 
X( I ,J) = 1.0 

3 CO^TINUF 

r PUT THF A MATRIX INTO ITS PLACE AT X ARRAY 

on 6 1 = 1 , m 
no a J = 1 , N 

A X( ( I +N ) • ( J + M ) ) =A ( I ,J) 

K3 I S THF COLUMN OF X WHERE F I R c T COLUMN n F THE A 
A NO 0 MATPICFS ARE 
on 21 J = X 3 , X 1 

RIG IS THF LARGEST ELFMFNT IN COLUMN 
FFAPCH INITIALISED TO Z c RO 
R I G=0 • 0 

AFTFR IC COLUMNS PF A IS D UT TO ITS NORMAL FORM RY 
PPW t RANSFPPMATIONS GO TO 22 
c IF ( ( J-M) .GT. IC) GO TP 2? 

K4=K2-X3+ J 
on 11 I =x a , x i 

SFARCH THF COLUMN OF A FOR L ARGF$ T FLFMENT fi^LOW K4 
MEMORIZE ROW AND COLUMN OP LARGEST FLPMENT 
IP(PIG.GP.DARS( X C 1 ,J) ) ) GO TO U 
R I G=DA P S ( X f I , J) ) 

TPOW=I 
JCOL= J 

11 CONTI NIJF 

TPf BIG. NF. 0.0) GO TO 1A 

IF ALL ELEMFNTS OF COLUMN ARE ZfRO PEMOVF THIS 
COLUMN SHIFT ALL COLUMNS AT RIGHT OF THIS ONE PLACE 
TO TH C LFFT AND PUT REMOVED COLUMN TO LAST COLUMN 
DO 13 11=1 , K 1 
TPMP=X( II , J) 

XK 1 =K 1— 1 

on l? j . i= j , x x i 

12 X( I I , JJ) = X( 1 I ,( JJ+l) ) 

X( I I ,X1 )=TFMP 

13 CONTI NIJF 
GO TO F 

1A IFf TROW.FO.XA) GO TO 16 

TP PIG IS NOT IN APPROPRIATE ROW CHANGE ROWS 
DP IE J 1 = 1 , K 1 
TFMP = X( XA, 11 ) 

X( K4,J1)=X( I ROW , J1 ) 

X( I ROW , J1 ) =TEMP 
IF CONTI NUF 

DIVIPF ROW OF BIG BY BIG 

16 TPMP = X< XA, J) 

no it j i = i . k l 

17 X(KA,J1)=X(KA,J1 )/TFMP 
DO 1® I =K 2 , K 1 
T P ( T « F 0 , K A ) GO TO 19 
IFf X< I ,J) .EO.O.O) GO TO 19 
G=Xf I , J) 

DO IP J1=1,X1 

X( I ,J1 )=X( I , J1 ) -G*X ( KA , J1 ) 

!F< OABSf X( I, J1 ) ) .LT.PD) X(I,J1)=0.0 
IP CONTINliF 
I'’ CONTI NUF 

21 CONTINUE 

AT t HIS PCI NT ALL RCWS OF TH c A BELOW ICTH SHOULD 
RF 7FR0 » IF RANK F QUA L TO COLUMN DIMENSION OF 
THF A, NORMAL FORM IS OBTAINED, OTHERWISE BY COLUMN 
tpANSFORMATI ONS FINISH the ALGOR I T H M 

22 X5=X3+IC 
K6=N+IC 
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L1=K3 

I f ( i r .fo. n> r,o to 41 

HP 33 l =*?,KC 
on 32 J=K5,K1 

IF f X( l , J ) .PO.0.0) r,0 TO 32 

0 = X ( L ♦ J ) 

on 31 I = 1 . K 1 

X( T , J) = X( f ,.I)-C*X( I ,L1 ) 

T F r 0A R S( X( I , J) ) . LT.DP) X(I , J)=0.0 



31 


CONTI NUF 








32 


CPM T IN)|F 










L 1 = 1 IM 








33 


CHMTI nijf 










PUT P 


TP 


ITS 


ARRAY 


I 


on ?3 1=1 


,M 








OP 23 J=1 


, W 






23 


T , J ) =X( 


( I 


T N ) , 


J) 




PUT o 


TO 


ITS 


ARRAY 




OP 2 L 1=1 


, N 








00 24 J=1 


, N 






24 


Q ( I , J ) = X ( 


I, 


( J+M) ) 




DF T ||P N 










EMP 









SUB POUT I NF mult 

CIIRCPCC ; 

SUBROUTINE MlItT FTMOS THF (IC,IM> NORTHWEST 
CCPNEP OF C=AB MATRIX M, ILTI PL ICAT T ON 

OFSopIPTirN OF PAPAMFT C RS: 

A — MULTIPLICAND MATRIX 
B — MULTIPLIER MATRIX 
C — RESULTANT matrix 

I. — ACTUAL COLUMN 01 Me NS T r N OF T H c A 
TC — ROW DIMENSION OF 0 E C T RFO NORTHWEST CORNER 
!M — COLUMN DIMENSION OF DESIRED N.W. CORNER 
A,B,L,IC»IM APE INPUTS, C IS OUTPUT 



SURROl’TINF MULTI L, A,B ,IC ,IM,C) 

IMc L IC I T R F A L*B ( A-H , D - 7 ) 

01 MENS I CM A (20,20) , B ( 2C , 20 ) , C ( 20 , 20 ) 
00 3 ! =1 » IC 

OO ? J= 1 ,IM 

RU M =C. C 
DO I K = 1,L 

1 «um = «;i)M+MI ,K)*B(K,J) 

2 C( I , J )=«UM 

3 CONTINUE 

RE T()P \t 
cKir 
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